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EXECUTIVE SUMMARY 


This report was sponsored under the Federal Aviation Administration's Aircraft 
Catastrophic Failure Prevention Research Program, Mission Need Statement No. 066- 
1 10. This report documents the NASA recommended Matched-Filter-Based One- 
Dimensional Search Method for the gust load analysis of nonlinear aircraft. The FAA 
Technical Program Monitor was Mr. Terry Barnes, ANM-105N, National Resource 
Specialist, Flight Loads and Aeroelasticity. FAA COTR was Mr. Thomas DeFiore, ACD- 
220, Principal Investigator, Flight Loads. 
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INTRODUCTION 


This manual and the accompanying code were developed as partial fulfillment of a 
NASA agreement with the FAA. The agreement had three main tasks. First, two 
NASA developed gust analysis methods were to be brought to the same level of maturity. 
These analysis methods were the Matched-Filter-Based (MFB) method and the 
Stochastic-Simulation-Based (SSB) method. Upon completion of the development work, 
the second task was to compare the methods and make a recommendation selecting which 
approach was best suited for nonlinear analyses. This work was completed and the 
results given in a presentation at the Gust Specialists Meeting in LaJolla, California on 
April 22, 1993. At the meeting, it was recommended by NASA that the MFB one- 
dimensional search was the method of choice for time-correlated gust loads analysis of 
aircraft with nonlinear systems. The third and final task was to develop a transportable 
computer program and accompanying documentation for using the recommended 
method. 

This manual describes a computer code called "MFB IDS" which performs the Matched- 
Filter-Based one-dimensional search. The MFB one-dimensional search is a 
deterministic method for obtaining maximized and time-correlated design loads for 
aircraft with nonlinearities. This method can, however, be applied to both linear and 
nonlinear aircraft. The paper summarizes the method, discusses the selection of gust 
intensity for the method, describes the FORTRAN code MFB IDS that performs the 
calculations, and presents numerical results for an example aircraft. 

BACKGROUND 

With the advent of aircraft that contain large numbers of nonlinearities in their flight 
control systems and/or gust load alleviation systems, existing methods for certifying 
aircraft for gust loads may not be adequate. For several years NASA Langley Research 
Center has conducted research in the area of time correlated gust loads and has published 
a number of papers on the subject (refs. 1-6). The initial research was restricted to 
mathematically linear systems (refs. 1-3). Recently, however, the focus of the research 
has been on defining methods that will compute design gust loads for an airplane with a 
nonlinear control system (refs. 4-6). To date, two such methods have been defined: one 
is based on matched filter theory; the other is based on stochastic simulation. 

The Matched-Filter-Based (MFB) method was developed first and was reported on in 
reference 4. The MFB method employs optimization to solve for its answers and this 
method comes in two varieties: the first uses a one-dimensional search procedure and the 
second a multi-dimensional search procedure. The first is significantly faster to run and, 
based on experience with a number of nonlinear models gives design loads only slightly 
lower in magnitude than the second. 
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The Stochastic-Simulation-Based (SSB) method has evolved over the past two years. 
References 5 and 6 describe the method. In reference 6 a comparison of the MFB and 
SSB methods was made. The results predicted by the two methods are strikingly similar 
and demonstrate that the key quantities from the MFB method (viz. critical gust profile, 
maximized load, and time-correlated load) are realizable in a stochastic analysis. 

Another significant finding in reference 6 was the relative computational costs of 
performing analyses using tf e MFB and SSB methods. For linear models, the MFB 
method is much more efficient than the SSB method. For nonlinear models, the MFB 
multi-dimensional search is much more expensive that the SSB method, while the MFB 
one-dimensional search requires less time that the SSB method. 

Since the MFB multi-dimensional search method as now implemented is prohibitively 
expensive, the options for methods that can practically be applied to nonlinear systems 
are the MFB one-dimensional search and the SSB methods. 

Based on the results in reference 6, the one-dimensional search is able to predict the 
maximized loads for nonlinear systems. In addition, it requires less computer resources 
than the SSB method. These factors show the utility of using the MFB one-dimensional 
search as a means of obtaining time-correlated gust loads for aircraft with nonlinear 
control systems. 

REVIEW OF THE MFB METHOD FOR LINEAR SYSTEMS 

The purpose of this section of the manual is to review the basic matched filter concepts. 
A detailed theoretical development of the MFB linear method can be found in reference 
2. The signal flow diagram in figure 1 outlines the implementation and illustrates the 
intermediate and final products of the method. 

Transfer-function representations of atmospheric turbulence and airplane loads are 
combined in series and represent the "known dynamics" boxes in the figure. A transfer- 
function representation of the von Karman spectrum in chosen for the gust filter. Load y 
is the load to be maximized. Loads zj through z n are the loads to be time correlated with 
load y. There are three major steps in the process: 

Step i The application of an impulse function of unit strength to the 

combined linear system, producing the impulse response of load y, 
h(t). 

Step ii The normalization of this impulse response by the square root of its 
energy and then reversing it in time. 

Step iii The application of this normalized reversed signal to the combined 
linear system, producing time histories of load y and time histories 
ofloads zj through z n . 

For simplicity of discussion these three steps will be referred to as the "MFB linear 
method." 
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The square root of the energy is defined by 



where h(t) is the impulse response of load y. For a stable system the numerical value of 
the quantity on the right side of equation (1) approaches a constant value as t 0 is 
increased. 

The selection of t 0 is based on the time required for the load impulse responses to damp 
out to near zero. Figure 1 shows that the impulse response of load y has essentially 
damped out at time t 0 . Thus, to is large enough for the response shown in figure 1 . Too 
large a t G value will, however, unduly increase the amount of computations required. The 
analyst must choose a value that provides accurate results while minimizing the use of 
computer resources. 

Within the time history of load y in step iii, the maximum value is y max . For linear 
systems, the theory guarantees that no other signal similarly normalized will produce a 
value of y larger than y, nax . This guarantee is a fundamental result of the MFB linear 
method. 

MFB ONE-DIMENSIONAL SEARCH DESCRIPTION 

The goal of Matched Filter Theory as applied to nonlinear systems is the same as that for 
linear systems: to find the maximized response time history, the maximum value of the 
response within that time history, and the time-correlated response time histories. 

Because the systems are not linear, the superposition principle of the MFB linear method 
no longer holds and the solutions for maximized loads cannot conveniently be obtained 
directly. The only practical means of finding the excitation waveform that maximizes 
ymax is a search procedure. The search is conducted systematically, subject to the 
constraint that the excitation waveform have a "unit" energy. 

Because superposition no longer holds, the magnitude and character of the responses are 
not necessarily proportional to the magnitude of the input. For the remainder of this 
paper two input magnitudes are important: k, the strength of the initial impulse; and o g , 
the design value of the gust intensity. (For the MFB linear method, the magnitude of 
both of these quantities was unity.) For nonlinear systems and a specific ag, the shape of 
the excitation waveform is a function of k , and, consequently, the quantity y max is also a 
function of this parameter. 

The one-dimensional search procedure performs a systematic variation of the quantity k 
to find the shape of an excitation waveform that maximizes y max . Figure 2 contains a 
signal flow diagram for this search procedure. Figure 2 is very similar to figure 1, but 
contains some subtle yet important differences that are indicated by the shaded boxes and 
by quotation marks. In figure 2 the initial impulse has a non-unity strength; the aircraft 
loads portion of the known dynamics box contains nonlinearities; and the shape of the 
excitation waveform and the value of y max are functions of the initial impulse strength. 

In addition, the "matched" excitation waveform and the "matched” load are shown in 
quotes because, for nonlinear systems, there is no guarantee that y max is a global 
maximum. 
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The application of the one-dimensional search procedure is as follows: 


Step 1 Select a specific design value for a g . 

Step 2 Select a set of values for k. 

Step 3 For each value of k, perform steps i through iii of the MFB linear 
method to obtain a set of "ymax" values and also the corresponding 
"matched" excitation waveforms. 

Step 4 From step 3 above, find the maximum value of y max and its 
corresponding "matched" excitation waveform. 

As seen in the first pass (upper half) of figure 2, the variation of the impulse strength (k) 
affects the MFB one-dimensional search analysis by changing the shape of the excitation 
waveform. For sufficiently low impulse strengths, the shape of the excitation waveform 
fo T nonlinear models will be invariant with k. While in this invariant region, the 
excitation waveforms will be similar to those that are obtained for linear models. For 
larger intensities the system nonlinearities are engaged and will cause the impulse 
responses and corresponding excitation waveforms to change shape. Consequently, they 
will differ from those obtained by linear systems. 

As seen in the second pass (lower half) of figure 2, the gust intensity affects the MFB 
one-dimensional search analysis by scaling the excitation waveform prior to being 
applied to the nonlinear model. Consequently, a low gust intensity should result in the 
nonlinear model behaving linearly. As gust intensity is increased beyond some threshold 
the nonlinear model response will begin to deviate from that of its linear counterpart. 

The effects of varying impulse strength and gust intensity will be discussed further in the 
numerical results section of the paper. 

SELECTION OF GUST INTENSITIES 

The purpose of this section of the paper is to present the reasoning behind the selection of 
the values of o g . Reference 6 describes the selection of gust intensities for the MFB and 
SSB methods. 

The following equation, from reference 7, expresses the "design value" of quantity y as 
defined in the design envelope criterion 

Vj.^n = KU a ( 2 ) 

where the quantity Ay is the RMS value of quantity y per unit RMS gust intensity, 
obtained from a conventional random process analysis of the airplane and U 0 is specified 
in the criterion. From reference 8 the quantity U 0 in equation (2) is shown to be the 
product of the gust RMS value and the design ratio of peak value of load to RMS value of 
load, and, therefore the quantity y design is interpreted as a peak value. 
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Reference 4 shows that, as a consequence of the normalization of the excitation 
waveform by the square root or its own energy and the use of unity gust intensity, the 
quantity y ma x from the MFB linear method is equal to the quantity X y from a 

conventional random process analysis, or 

y m J0 g =1 ) = \ 0) 

In equation (3) ymax is interpreted as an RMS value, not a peak value. Substituting 
equation (3) into equation (2), yuesign is now 


y design ~ y»iax ( ® g 1)^0 

If, in performing the MFB linear method, U 0 is used for the gust intensity then the 
quantity y max is equal to 


y'liiax ( ® g — ) — Wo ( 5 ) 

The right hand sides of equations (2) and (5) are seen to be equal, therefore 

y design ~ V max ^ g ~ 

Two options for the value of cs g have been offered: = 1, for which y design is defined by 

equation (4); and a g = U a , for which y<jesign is define! by equation (6). When analyzing a 
linear system the choice of c g is irrelevant because the same value of y design will be 
obtained in either case. However, when nonlinearities are introduced into aircraft control 
systems, loads are not simply proportional to gust intensity. Consequently, o g should be 
set to U c in the MFB nonlinear calculations, or 

® k MFB = U a ( 7 ) 

and the resulting "y max " values from the method should be interpreted as y design- 


DESCRIPTION OF MFB1DS 

MFB IDS is a FORTRAN 77 program which performs the Matched-Filter-Based one- 
dimensional search. The MFB IDS solution procedure follows what was outlined in 
figure 2. The code must be run once for each combination of gust intensity and output 
quantity for which a maximized value is desired. 

Figure 3 shows the solution procedure used by MFB IDS. The program first reads the 
input data then generates the impulse responses for each k value by calling the simulation 
subroutine. The simulation subroutine uses a public domain ordinary differential 
equation solver to generate output time histories using a user defined subroutine 
containing the aircraft equations of motion. Next, the excitation waveforms are generated 
by reversing the impulse responses in time and normalizing them by the square root their 
respective energies. Then, the simulation subroutine is again called for each of the 
normalized excitation waveforms to obtain the "maximized" load response time histories. 
Finally, the pertinent output quantities are written to computer files. 
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This section of the manual describes the main parts of MFB1DS. 

Required Files 

Six files are required to run MFB IDS. 

MFB1DS.F 
MFB IDS. INC 

MODEL. F 

MODEL.INP 
LSODE.F 
INTUTILS.F 

The two public domain files LSODE.F and INTUTILS.F are well documented in their 
respective source codes and will not be discussed in detail in this manual. These files 
were selected because of the fact that they are in the public domain. The user is 
encouraged to substitute more efficient ordinary differential equation solvers if available. 

MFB1DS.F Subroutines 

Subroutine DATAIN 
Subroutine SIMULATE 

Subroutine SAVMATFORM 

Subroutine MATSAV 

Subroutine SAVASCIIFORM 

MODEL.F 

MODEL.F is a user supplied file containing the subroutine EQSMOT, the equations for 
the gust filter in series with the aircraft equations of motion. This subroutine uses 
variables x (vector of states) and u (excitation) to obtain xd (derivative of x) and y (vector 
of outputs). 


Reads the input parameters from the input file 

Computes output time histories using LSODE.F 
and MODEL.F 

Saves output quantities to a MATRIX; _ readable 
file 

Used by Subroutine SAVMATFORM to save data 
in MATRIX x format 

Saves output quantities in a more easily readable 
ascii file 


Main program, see Appendix A 

File containing the common blocks used by 
MFB IDS. F, see Appendix B 

User supplied file containing the 
subroutine EQSMOT, the aircraft equations of 
motion, see Appendix C 

User supplied file containing the input 
quantities required to run the code, see Appendix D 

Public domain ordinary differential equation 
solver 

Public domain subroutines used by LSODE.F 
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The following lines of code must be present at the top of the file MODEL.F: 

subroutine EQSMOT(neq,t,x,xd) 
parameter(maxout=20) 
double precision x(*), xd(*), y(maxout) 
double precision t, tstart, tend, uO, ul 
common /eqsmotcom/y,tstart,tend,uO,u 1 
c 

u=(t-tstart) *(u 1 - u())/(tend- tstart)+u( ) 

In spite of the fact that the user specifies a desired time step, the ordinary differential 
equation solver LSODE calls the EQSMOT subroutine at values of time (t) between time 
steps. Consequently, the equation on the last line of the above code is included to 
provide a linear interpolation for values of u between time steps. 

A linear system can be used to demonstrate how x, xd, y and u are related. For a linear 
system only, the equations of motion can be written in state space form: 

{xd} = [A]{x} + [B]u 
{y} = [C]{x} + [D]u 

While equation (8) is linear, the relationship between x, u and xd is, in general, a 
nonlinear one. 

Description Of Impulse Input 

Since the impulse generating procedure is one of the key components of the program and 
there are several methods that can be used to generate impulse functions, a few words 
describing the impulse function procedure implemented in the program are warranted. 
The procedure chosen to generate the impulse function in this code is straightforeword 
and can be found in main program listing in Appendix A. In general, the impulse input 
function, as seen by the differential equation solver, is a ramp up from zero at the first 
time step to a value of k/(2*deltat) at the second time step, a constant value of 
k/(2*deltat) between time steps two and three, and a ramp down to zero between time 
steps three and four. 

♦ Description Of Output 

Key information is written to standard out (unit 6) as shown in Appendix E. This 

* information allows the user to monitor the progress of the program. Two additional 
output files can be created which contain time history and correlated output information. 
Both of these output files are ascii files: one is in a MATRIXx readable form and the 
other is in a more easily read form. The variables contained in the MATRIXx readable 
file are discussed below. The quantities found in the other output file are self explanatory 
and will not be discussed here. 

The following scalar quantities are written to the MATRIXx file: 
sigmag Gust intensity, 

noutmx The output quantity to be maximized. 


7 


deltat 


The lime step. 

tmaximp The length of the impulse responses. 

The following arrays can also written to the MATRIX x file: 

allkvals This vector has length nkvals, and stores all the k values used in 
the analysis. 

maxout The maximum output values for each k value are stored in each 

column. Thus, this array will he an nout by nkvals array. 

impres# An array of this form is created for each output quantity. For 

instance, for output 1 an array named "impres 1" will be created and 
for output 5 an array named "impres5" will be created. Each 
column of these arrays are impulse response time histories for one 
of the k values. Thus, these are ntsteps by nkvalues arrays. 

wavef# An array of this form is created for the output to be maximized. 

For instance, if output 6 is to be maximized, then an array named 
"wavef6" will be created. Each column of this array is an 
excitation waveform matched to output # for a k value. Thus, this 
is an ntsteps by nkvalues array. 

exresp# An array of this form is created for each output quantity. For 
instance, for output 1 an array named "exresp 1” will be created. 
Each column of this array is an excitation waveform response for a 
k value. Thus, these are (2*ntsteps+l) by nkvalues arrays. 

NUMERICAL EXAMPLE 

This section describes the step by step process of obtaining numerical results at a gust 
intensity of 1,530 in/sec. 

Step 1 - Create a subroutine containing the gust filter in series with the aircraft 
equations of motion. 

The nonlinear simulation model of the ARW-2 drone aircraft equipped with a nonlinear 
control system was used in this example. This model was connected in series with a 
transfer-function representation of atmospheric turbulence. The transfer function used in 
this example was (ref. 8) 

*g_ _ a PL [1 + 2.618(L/V)s][l + 0.1298(LJV)s] 

r] ° 8 V nV [1 + 2. 083(L/V)s ][! + 0. 823(L/V)s ][1 + 0. 0898(L/V)s ] 

which approximates the square root of the von Karmon power spectral density function. 
The quantity o g is the intensity of the gust or standard deviation — which, assuming zero 
mean, is also equal to the root-mean-square, or RMS, value - of gust velocity. 

Figure 4 shows a block diagram of the ARW-2 simulation model and includes the 
aeroelastic plant, a gust load alleviation (GLA) control law, and nonlinear control 
elements. The aeroelastic plant is a linear, s-plane aeroelastic half-model consisting of 
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two longitudinal rigid-body modes and three symmetric flexible modes. Unsteady 
aerodynamics were obtained using the doublet lattice method (ref. 9). The model also 
includes the dynamics of the control surface actuators. The two-input/two-output GLA 
control law was obtained using a Linear Quadratic Gaussian design approach with the 
intent of reducing wing root bending moment (ref. 10). The nonlinear elements impose 
deflection limits of ±1" on the elevators and 0“ to +1" on the ailerons to simulate spoilers 
The model contains 32 states, and the analysis conditions are at a Mach number of 0.86 
and an altitude of 24,000 feet. 

The resulting combined system has a single input (white noise) and many outputs 
including the gust velocity (y(14)). The output quantities for this model are as follows: 


y(D-y(3) 

= No physical interpretation 

y(4) 

= Tip acceleration 

y(5) 

= Fuselage acceleration 

y(6) 

= Wing root bending moment (WRBM) 

y(7) 

= Wing root shear (WRS) 

y(8) 

= Wing outboard bending moment (WOBBM) 

y(9) 

= Wing outboard torsion moment (WOBTM) 

y(io) 

= Elevator deflection 

y(ii) 

= Elevator rate 

y(i2) 

= Aileron deflection 

y(i3) 

= Aileron rate 

y(14) 

= Gust velocity 


Appendix C contains the FORTRAN version of the analytical model. This model was 
created using MATRIX x SYSTEM BUILD (ref. 1 1) and converted to FORTRAN using 
the HYPERCODE (ref. 12) package. As shown in Appendix C, the appropriate lines of 
code were added to the model created by the HYPERCODE package to make it 
compatible with the differential equation solving scheme used in MFB1DS. 

Step 2 - Create an input file for the model. 

The length of the state vector and number of model output quantities are contained in the 
input file MODEL.INP and must be made compatible with the analytical model created 
in step 1. Other input quantities like the gust intensity, identifying the output quantity to 
be maximized, and the range of k values must also be selected. (Trial and error will be 
required to find the range of k values where the maximum load value is obtained.) In 
addition to these parameters, the length of time for the impulse responses and the time 
step must be chosen. The length of the impulse responses should be only long enough to 
allow the output quantity to damp out to a relatively small value as shown in the example 
in figure 4 (a), while the time step should be chosen in accordance with the model 
response characteristics. 

The input file used in the ARW-2 analysis is shown in Appendix D. 

Step 3 - Modify the parameters. 

Appendix B contains a listing of MFB 1DS.INC and a description of the parameters and 
variables contained there. The parameters found in the file MFB IDS. INC and MODEL.F 
must be made large enough to accommodate the analytical model (step 1) and the input 
file (step 2). Computer storage can be minimized by making parameters "maxstates" and 
"maxout" equal to the minimum values required for a given analytical model. 
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Step 4 - Compile Ihe code. 

The following statement compiles the code on a SUN workstation. 
(77 MFB I DS.F MODEL. F LSODF..F INTUTILS.F 


Slop 5 - Execute the code. 

Appendix D shows (he initial input file used in the analysis of the ARW-2 model. The 
initial value of gust intensity was 1.530 in/sec (1.5X85 ft/scc). This initial value of gust 
intensity was chosen so that the nonlinearities would he invoked in this example. Nine k 
values were used ranging from 10 to I5,0(X). The output quantity to be maximized was 
y(6). wing root bending moment. WRBM. The standard output from this run is shown in 
Appendix E. 

The following statement will run the code on a UNIX system 
a.out < MODEL. INI’ 

Step 6 - Cxnniine output. 

Of the nine k values used in the analysis, all the intermediate results for three of the nine 
k values arc shown in figure 5. Beginning with plot (a), the impulse responses for the 
three k values are shown. Plot (b) gives the corresponding excitation waveforms obtained 
from these impulse responses. Plot (c) shows the critical gust profiles and plot (d) shows 
the maximized load responses created by the gust input. Plot (e) depicts maximum 
WRBM as a function of initial impulse strength: the circles are actual numerical results; 
the solid line is a faired line. For this particular example, a value of k of about 2.410 
creates an excitation waveform and critical gust profile that yields the largest maximized 
value of WRBM. 296.994 in-lbs. 

Step 7 - Vary parameters. 

The user may wish to refine the results shown in figure 5. With this in mind, a new range 
of k values was selected and the case was rerun. Figure 6 shows the results of using the 
new set of k values distributed between 400 and 6,000 producing a much smoother curve. 
The maximum load value obtained using the refined k distribution was 296,804 in-lbs at a 
k value of approximately 2,1 73. For this particular airplane and set of flight conditions a 
larger load value was not obtained with the new range of k values. 

EXAMINATION OF RESULTS 

This section has been included in the manual to provide the user with insight into how the 
results Irom Ml BIDS should be interpreted. Ihe results presented in the previous 
section will be shown along \^ith results obtained using two additional gust intensities. 

To obtain results at different gust intensities, the sigma value of the input File shown in 
Appendix D was changed and MFB IDS was rerun. In addition, answers were obtained 
for each of the three gust intensities using a linearized version of the model. For 
comparison these linear answers will be plotted with the results from the nonlinear 
model. Note that when using a linear model, only one k value needs to be used because 
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the answer will not be a function of k. Consequently, the maximum load value plotted 
versus k is a horizontal line for a linear model. 

The results for gust intensity values of 1,020 in/sec, 1,530 in/sec and 2,040 in/sec are 
shown in figure 7. Before proceeding, it is important to point out that results are problem 
dependent, and while these results are typical of all the aircraft examined by NASA, 
different aircraft with different nonlinearities may exhibit different trends. With this in 
mind, two important points will be made concerning the general trends exhibited in figure 
7. First, regardless of the gust intensity value, the maximum load is constant for small 
values of k. Consequently, when searching for the maximum load, the range of interest 
for the k values can be limited on the low end. Second, the maximum load value 
decreases toward some relatively small value as k is made very large. Thus, the range of 
interest for k values can also be limited on the high end. What happens in between these 
two extremes depends on the specific gust intensity. 

In figure 7 the maximum load values are constant for k less that 400. In plot 7(a) a 
significantly larger load value was not found when k was increased beyond 400, while for 
plots 7(b) and 7(c) a larger value was found. This indicates that the character of the 
results (i.e., the shape of the maximum load versus k plot) and the specific k value that 
produces the maximum load value are functions of gust intensity. But, the range of k 
values where the maximum will be found is generally limited, and in this case the range 
is roughly between 400 and 10,000. It should also be noted that the difference between 
answers from linear and nonlinear models is also a function of the gust intensity. Larger 
gust intensity values generally result in larger differences between answers from linear 
and nonlinear models. Here, this difference is 2% for the lowest gust intensity and 18% 
for the largest gust intensity. 


CONCLUDING REMARKS 

This manual has reviewed the theory behind the Matched-Filter-Based one-dimensional 
search procedure. The code that performs this procedure has been discussed and example 
numerical results were presented and interpreted. The code, MFB1DS, is available in a 
self contained form. It has all the required equation solvers and files necessary to run the 
example problem. The user is, however, encouraged to modify the existing code by 
inserting more efficient routines if available. 
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Figure 1 . MFB linear method signal flow diagram. 
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Figure 2. Nonlinear MFB signal flow diagram for the one-dimensional search. 
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Figure 3. MFB IDS solution procedure. 
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Figure 4. Nonlinear block diagram for the ARW-2 drone aircraft. 
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Appendix A - MFB1DS.F source code listing 
program MFB1DS 

c************************************************************ 

c* 

c* MFB1DS.F: 
c* 
c* 

C* 

c* written by: 
c* 
c* 
c* 
c* 
c* 

c* compile: 
c* 
c* 

£************************************************************ 
include MFBIDS.INC' 
character* 80 ctitle 
character* 80 cease 
common /case/ccase 
common /title/ctitle 
dimension yth(maxout,2*maxtsteps) 

c **read input parameters 
call DATAIN 

write(6,*) ctitle 
write(6,*) cease 

c **calculate number of time steps for impulse response 
ntsteps=int(tmaximp/deltat + 0.001) + 1 
write(,6,*) ’ ’ 

write(6,*) ’ length of simulation for impulse responses = 

+ tmaximp 

write(6,*) ’ number of time steps for impulse responses = 

+ ntsteps 
write(6,*) ' ' 

write(6,*) ' length of simulation for excitation responses = 

+ 2*tmaximp 

write(6,*) ’ number of time steps for excitation responses = 

+ 2*ntsteps- 1 

write(6,*) ' ' 

write(6,*) ' gust intensity = sigma 
write(6,*) ' ’ 

c **create array of zeros for inpulse input 
do 100 i=l,(2*ntsteps-l) 

100 aimpt(i)=0.0 

c * Calculate k values 
do 200 k=l,nkvals 


a matched- filter-based method of obtaining * 

maximized and time-correlated gust loads for * 

linear and nonlinear aircraft * 

* 

Robert C. Scott * 

NASA Langley Research Center MS/340 * 

Hampton, VA 23665 * 

* 

804-864-2838 * 

* 

177 MFB1DS.F MODEL.F LSODE.F INTUTILS.F * 

* 

where MODEL.F is the a/c model * 


A- 1 


if (nkvals.eq.l) then 
allkvals(l)=akmin 
else 

allkvals(k)= 10.0**( (k- 1 )*(log 10(akmax)-log lO(akmin))/ 

+ float(nkvals-l)+logl()(akmin) ) 

endif 

200 continue 

c **obtain impulse response for each k value 

write(6,*) 'calculating impulse responses for each k value' 
do 300 k=l,nkvals 
ak=allkvals(k) 

c **assign impulse strength k to array aimpt 
aimpt(2)=ak/deltat/2 
aimpt(3)=ak/deitat/2 

c **call subroutine SIMULATE to obtain impulse response 
call SIMULATE(yth,aimpt,ntsteps) 
do 300 iout=l,nout 
do 300 j=l,ntsteps 

300 aimp_res(j,iout,k)=yth(iout,j) 

c ’•‘♦calculate normalized excitation waveforms 

write(6,*) 'calculating normalized excitation waveforms' 
do 400 k=l,nkvals 
i=noutmx 
energy=0.0 

energy=aimp_res(l,i,k)**2 + aimpjres(ntsteps,i,k)**2 
do 401 j=2,ntsteps-l 

401 energy = energy+2 + (aimp_res(j,i,k))**2 

energy = sqrt( tmaximp+energy + ,5/float(ntsteps)/3. 14159 ) 
write(6,*) ' ', k, ' k = ', allkvals(k), 

+ ' sqrt(energy) of output', noutmx, ' = ', energy 

do 402 j=ntsteps, 1 1 

402 wave(ntsteps- j+ 1 ,k)=sigma*aimp_res(j,i,k)/energy 
400 continue 

c ♦'•'obtain maximized responses 

c Note that here the length of the simulation is 2*tmaximp and the 
c number of time steps for maximized responses is 2*ntsteps-l 

write(6,*) 'calculating maximized responses' 
do 500 k=l,nkvals 

do 501 j=l,ntsteps 

501 aimpt(j)=wave(j,k) 

do 502 j=(ntsteps+l), (2*ntsteps-l) 

502 aimpt(j)=0.0 

call SIMULATE(yth,aimpt, (2*ntsteps-l) ) 
do 503 j=l,(2*ntsteps-l) 
do 503 i=l,nout 

503 aexc_res(j,i,k)=yth(i,j) 


A- 2 


write(6,*) ' k,' k = \allkvals(k), 

+ ' maximum value of output 

+ noutmx, ' = aexc_res(ntsteps,noutmx,k) 

500 continue 

c **save data 

write(6,*) 'saving data' 

if (iouttype.eq.l.or.iouttype.eq.3) call SAVASCIIFORM 
if (iouttype.eq.2.or.iouttype.eq.3) call SAVMATFORM 

stop 

end 


cl SIMULATE: subroutine to obtain a time history for the a/c model I 

c 

subroutine SIMULATE(yth,u,nsteps) 
include 'MFBIDS.INC' 

dimension yth(maxout,2*maxtsteps), u(2*maxtsteps) 

double precision y(maxout) 

double precision atol, rtol, xstate(maxstates) 

double precision rwork(22*l()*maxstates+(2*l+l)’ , ‘maxstates) 

double precision tstart.tend, u(),ul 

integer neq,itol,itask,istate,iopt,lrw,liw,mf 

integer iwork(20+maxstates) 

common /eqsmotcom/y,tstart,tend,u(),u 1 

c **for guidance in selecting sizes for arrays rwork and iwork 
c see source code LSODE.F 

external EQSMOT 

c **Setup input for LSODE 

c see LSODE.F for explanation of these parameters 

c 

neq=nstates 
c tstart=starting time 
c tend=ending time 
itol=l 

rtol=1.0E-06 

atol=1.0E-10 

itask=l 

istate=l 

iopt=0 

c rwork=real work space 

lrw=22* 10*maxstates+(2* 1 + 1 )*maxstates 
c iwork=imiginary work space 
liw=20+maxstates 
c iaceqs 
mf=23 

do 1 i=l,maxstates 
1 xstate(i)=0.() 


n o 


do 1 1 i=l,maxout 
11 yth(i,l)=0.() 

tstart=0.0 

tend=0 

do 2 i=2,nsteps 
tend=tend+deltat 
uO=u(i-l) 
ul=u(i) 

call LSODE(EQSMOT,neq,xstate,tstart,tend,itol,rtol,atol, 
+ itask,istate,iopt,rwork,lrw,iwork,liw,jaceqs,mf) 
c **If istate is not equal to 2 then an error has occurred 
c and the following will be printed 

if (istate.ne.2) write(6,*) 'istate = istate, 

+ 'tend = ', tend 

tstart=tend 
do 3 j=l,nout 
3 yth(j,i)=y(j) 

2 continue 

return 

end 


I DAT AIN: subroutine to read data from file MFB1DS.INP 1 

c 

c 

subroutine DATAIN 
include 'MFB1DS.INC 
character* 80 ctitle 
character* 80 cease 
character* 15 edata, cmatrixx 
common /case/ccase 
common /title/ctitle 
common /fnames/edata, cmatrixx 
character cdummy*80 

c **read model title 

read(unit=5,fmt='(a80)') cdummy 
read(unit=5,fmt=’(a80)’) ctitle 

c **read case title 

read(unit=5,fmt=’(a80)') cdummy 
read(unit=5,fmt='(a80)') cease 

c **read nstates, nout, noutmx 

read(unit=5,fmt='(a80)’) cdummy 
read(unit=5,fmt=*) nstates.nout, noutmx 

c **read tmaximp, deltat, nsubintvl 
read(unit=5,fmt='(a80)') cdummy 
read(unit=5,fmt=*) tmaximp, deltat 
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c **read sigma, akmin, akmax, nkvals 
read(unit=5,fmt='(a80)') cdummy 
read(unit=5,fmt=*) sigma, akmin, akmax, nkvals 

c ’•"•‘read iouttype, impres, ixcres, iwave 
read(unit=5,fmt='(a8())’) cdummy 
read(unit=5,fmt=*) iouttype, impres, iexcres, iwave 

c **read data file name 

read(unit=5,fmt=’(a80)') cdata 
read(unit=5,fmt='(al2)') cdata 

4 

c **read matrixx file name 

read(unit=5,fmt=’(a80)') cdummy 
read(unit=5,fmt='(al2)') cmatrixx 

close(unit=5) 

return 

end 

c - - - - 

cl SAVMATFORM: subroutine to write output to a MATRIXx readable rile I 

c - - 

subroutine SAVMATFORM 
include 'MFB1DS.1NC 
dimension rtemp(2*maxtsteps,maxstates) 
character* 10 varname, cnout 
character* 15 cdata, cmatrixx 
common /fnames/cdata, cmatrixx 

open(unit= 1 ,file=cmatrixx) 
rewind(l) 

c ** write gust intensity 
rtemp(l,l)=sigma 

call MATSAV(l,'sigmag', 2*maxtsteps, 1, 

, + 1, 0, rtemp, rtemp, '(Ip2e24. 15)') 

c **write k values for each case 
» do 1 i=l, nkvals 

1 rtemp(i,l) = allkvals(i) 

call MATSAV(l,'kvals', 2*maxtsteps, nkvals, 

+ 1,0, rtemp, rtemp, '(Ip2e24.15)') 

c **write which output quantity was "maximized" 
rtemp(l,l)=noutmx 

call MATSAV(l,'noutmx', 2*maxtsteps, 1, 

+ 1,0, rtemp, rtemp, '( 1 p2e24. 1 5)') 

c ** write time step 
rtemp(l,l)=deltat 

call MATSAV(l,'deltat', 2*maxtsteps, 1, 

+ 1, 0, rtemp, rtemp, '(Ip2e24. 15)') 
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c ** write tmaximp 
rtemp( 1 , l)=tmaximp 

call MATSAV( l, 'tmaximp', 2*maxtsteps, 1, 

+ 1,0, rtemp, rtemp, '(Ip2e24.15)’) 

c **write maximum output values 
do 3 i=l,nout 
do 3 k=l,nkvals 

3 rtemp(k,i)=aexc_res(ntsteps,i,k) 

call MATSAV(l,'maxout’, 2*maxtsteps, nkvals, 

+ nout, 0, rtemp, rtemp, '(Ip2e24. 15)') 

c ** write impulse responses 
if (impres.eq.l) then 
do 4 i=l,nout 

open(unit=50, status='SCRATCH') 
write(unit=50,fmt= '(i3)') i 
rewind(50) 
ndigits=l 

if (i.gt.9) ndigits=2 
if (i.gt.99) ndigits=3 
cnout=' ' 

read(unit=5(),fmt='(a3)') cnout(l:ndigits) 
close(50) 

vamame='impres'//cnout 
do 5 k=l, nkvals 
do 5 nt=l,ntsteps 

5 rtemp(nt,k)=aimp_res(nt,i,k) 

call MATSAV(1, vamame, 2*maxtsteps, ntsteps, 
+ nkvals, 0, rtemp, rtemp, '(Ip2e24.15)') 

4 continue 
endif 

c **write excitation responses 
if (iexcres.eq.l) then 
do 6 i=l,nout 

open(unit=50, status='SCRATCH') 
write(unit=50,fmt= '(i3)') i 
rewind(50) 
ndigits=l 

if (i.gt.9) ndigits=2 
if (i.gt.99) ndigits=3 
cnout=' ' 

read(unit=50,fmt='(a3)') cnout( 1 :ndigits) 
close(50) 

varname='exresp'//cnout 

do 7 k=l, nkvals 
do 7 nt=l,(2*ntsteps-l) 

7 rtemp(nt,k)=aexc_res(nt,i,k) 
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call MATSAV(1, vamame, 2*maxtsteps, (2*ntsteps-l), 
+ nkvals, 0, rtemp, rtemp, '( Ip2e24. 15)’) 

6 continue 
endif 

c ** write excitation waveforms 
if (iwave.eq.l) then 

open(unit=50, status='SCRATCH') 
write(unit=50,fmt= ’(i3)') noutmx 
rewind(50) 
ndigits=l 

if (noutmx.gt.9) ndigits=2 
if (noutmx. gt.99) ndigits=3 
cnout=' ' 

read(unit=50,fmt='(a3)') cnout(l:ndigits) 
close(50) 

vamame='wavef//cnout 
do 9 k=l, nkvals 
do 9 nt=l,ntsteps 
9 rtemp(nt,k)=wave(nt,k) 

call MATSAV(1, vamame, 2*maxtsteps, ntsteps, 

+ nkvals, 0, rtemp, rtemp, '( 1 p2e24. 1 5)') 

endif 

close(l) 

return 

end 


SAVASCIIFORM: subroutine to write output to an scii file I 


subroutine SAVASCIIFORM 
include ’MFB1DS.INC 
, character* 80 ctitle 

character* 80 cease 
character* 15 edata, cmatrixx 
» common /case/ccase 

common /title/ctitle 
common /fnames/edata, cmatrixx 

open(unit=9, file=cdata) 
rewind(9) 

c **write general information 
write(9,'(a80)') ctitle 
write(9,'(a80)') cease 
write(9,*) ' sigma = ', sigma 
write(9,*) ' matched output quantity = ', noutmx 
write(9,*) ’ tmaximp = ', tmaximp 
write(9,*) ' deltat = ', deltat 
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c **write maximum output values for each case 
write(9,*) ' ' 

write(9 *) '***************************************************' 
write(9,*) ‘MAXIMIZED AND TIME CORRELATED MAX OUTPUT 
QUANTITIES’ 

write(9,’ l ‘) '***************************************************' 
do 1 i=l,nout 

write(9,*) ’ output quantity = i 
write(9,*) ’ k value maximum output value' 
do 1 k=l,nkvals 
write(9,*) allkvals(k), ' 

+ aexc_res(ntsteps,i,k) 

1 continue 

c **write impulse responses 

if (impres.eq. 1) then 
do 2 k=l,nkvals 
write(9,*) ’ ’ 

write(9,*) '***********************************************' 
write(9l*) ’IMPULSE RESPONSES’ 

write(9 *) '***********************************************' 
write(9,*) ' kvalue = allkvals(k), 

+ ' maximized output quantity = noutmx 

do 2 iout=l,nout 

write(9,*) 'output quantity = iout 
do 2 j=l,ntsteps 
write(9,*) aimp_res(j,iout,k) 

2 continue 
endif 

c **write excitation waveforms 
if (iwave.eq.l) then 
do 3 k=l,nkvals 
write(9,*) ' ’ 

write(9 *) '***********************************************' 
write(9*) 'EXCITATION WAVEFORMS' 

write(9,*) ’***********************************************' 
write(9,*) ' kvalue = allkvals(k), 

+ ' output quantity = noutmx 

do 3 j=l,ntsteps 
write(9,*) wave(j,k) 

3 continue 
endif 


**write excitation responses 
if (iexcres.eq.I) then 
do 4 k=l,nkvals 
write(9,*) ' ' 

write(9 *) '***********************************************' 
write(9i*) 'EXCITATION WAVEFORM RESPONSES' 
write(9 *) '***********************************************' 
write(9,*) ' kvalue = allkvals(k), 
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+ ’ maximized output quantity = noutmx 

do 4 iout=l,nout 

write(9,*) 'output quantity = iout 
do 4 j=l,2*ntsteps-l 
write(9,*) aexc_res(j,iout,k) 

4 continue 
endif 


close(9) 

return 

end 


I MATSAV: write variables to a file in matrixx format I 

c 

subroutine MATSAV ( lunit, name, nr, m, n, img, 

+ xreal, ximag, formt ) 

c 

c 

c 

c MATSAV writes a matrix to a file in a format suitable for the 
c matrixx load operation, 
c 

c — 

c param. type on input- on output- 

c 

c 

c lunit integer fortran logical unit number. unchanged, 
c 

c name character^*) name of the matrix, one al- unchanged. 


c 

c 

p 

(maximum 
length 10) 

phabetic followed by up to 9 
alphanumeric characters. 


V 

c nr 

c 

c 

c 

c 

c 

integer 

row-dimension in the 
defining dimension or type 
statement in the calling 
program, nr must be greater 
than or equal to m. 

unchanged. 

c m 

c 

integer 

number of rows of the matrix 

unchanged. 

c n 
c 

c * 

integer 

number of columns of the 
matrix. 

unchanged. 

v> 

c img 

c 

c 

integer 

if img = 0, the imaginary 
part (ximag) is assumed to 
be zero and is not saved. 

unchanged. 

c 

c xreal 
c 

double 

precision 

real part of the matrix to 
be saved. 

unchanged. 

c ximag double 

imaginary part of the matrix 

unchanged. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


precision to be saved. 

formt character^*) string containing the for- unchanged, 
(maximum tran format to be used for 
length 20) writing the elements of the 
matrix. 


example: the following fortran program generates an elementary 

matrix in x and writes it to fortran unit 1 . assume 
that unit 1 has been preallocated as file (data set) 
test. 


dimension x(20,3), dummy 
do 200 j= 1,3 
do 100 i= 1,10 
x(i,j)=0.0d() 

100 continue 
x(j,j)=1.0d0 
200 continue 

call MATSAV( 1, 'amatrix', 20, 10, 3, 0, 
$ x, dummy, '(Ip2e24.15)' ) 

stop 
end 


after this program runs, invoke matrixx and type: 

<> load 'test' 

this will put x on the stack as stack-variable-name amatrix. 


integer lunit, m, n, nr, img 
character*^) name, formt 
dimension xreal(nr, 1 ), ximag(nr, 1 ) 

character nam * 1 0, form * 20 


write header record. 


nam=name 

form=formt 

write(lunit,'(al(),3i5,a20)') nam, m,n,img, form 


write real-part of the matrix. 


write(lunit,form) ((xreal(i,j),i= 1 ,m),j= 1 ,n) 


write imaginary-part if nonzero. 




t 
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if(img.ne.O) write(lunit,t'orm) ((ximag(i,j),i=i,m),j=l,n) 

return 

end 




» 
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Appendix B - MFB1DS.INC source code and description of parameters 

parameter(maxtsteps=5(X)l, maxstates=40, maxout=2(), maxkvals=20) 
common /modelinfo/nstates, nout 
common /siminfo/tmaximp, deltat, ntsteps, nsuhintvl 
common /case_info/noutmx 

common /excitation_info/sigma, akmin, akmax, nkvals, 

+ alikvals(maxkvals) 

common /t_hist/aimp_res(maxtsteps,maxout,maxkvals), 

+ aexc_res(2*maxtsteps,maxout,maxkvals), 

+ wave(maxtsteps.maxkvals), 

+ aimpt(2*maxtsteps) 

common /output/iouttype, impres, iexcres, iwave 


The following is a description of the parameters: 

maxtsteps This parameter is the maximum number of time steps for the 

impulse responses, (maxtsteps > tmaximp/deltat +1) 

maxstates This parameter is the maximum number of states in the equations 

of motion, (maxstates > nstates) 

maxout This the maximum number of output quantities required by the 

aircraft model, (maxout > nout) 

maxkvals This the maximum number of kvalues which can be run. 

To change these parameters edit MFB 1DS.INC. In addition, maxout is also found in 
MODEL.F and must be equal that the value in MFB IDS. INC. Computer storage can be 
minimized by making maxstates and maxout equal to the minimum values required for a 
given analytical model. 
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Appendix C - MODEL. F (ARW-2) source code listing 


subroutine EQSMOT(neq,t,x,xd) 
parameter(maxout=20) 
double precision x(*), xd(*), y(maxout) 
double precision t.tstart.tend, u(), u 1 
common /eqsmotcom/y,tstart.tend,u(),ul 

u=(t- tstart) * (u 1 - u())/(tend- tstart) +u0 


c - - state-space system 

} c -- {nlarw2.controller.l } 

c -- ss c 

y(l) = -5.244970395787696D-4*x( 1 ) - 2.747(X)469140994399D-4*x(2) 

• + - 1.42946149773359D-3*x(3) + 5.56258182659891597D-6*x(4) 

y(2) = -3.287044785655 14 102D-3*x(l) - 6.743878267402364D-5*x(2) 
+ - 9.87258 135859925306D-4*x(3) + 

+ 9.90556693932225 104D-6*x(4) 

c 1 - u bounded limit 

c — {nlarw2.aileron limiter.2} 

y(3) = min( 0.0D0, max( -0.0174499999999999999D0, y(2) ) ) 
c y(3) = min( 0.0 1 74499999999999999DO, max( 
c + -0.0 1 74499999999999999D0, y(2) ) ) 

c state-space system 

c - {nlarw2.arw2 moD.3} 

Q Q 

y(4) = -2846.47376 1074609D0*x(5) - 2699.5767200 1826D0*x(6) + 

+ 6978.2941 1601534196D0*x(7) + 2. 14350293 155(X)8902D5*x(8) - 

+ 22.328976390071 260 lD0*x(9) - 3.3990629437 14242D0*x( 10) - 

+ 9.13144768834393994D0*x(l 1) + 5.05173942454337099D0*x(12) 

+ + 75.7538 109829206405D0*x( 13) - 

+ 0.95882 175O950846804D0*x(14) - 

+ 0.86024628749677 13()3DO+x( 15) - 78.95852907402 1 5004D0 + x( 1 6) 

+ + 43.922 1734285092698D0*x( 17) + 

+ 8.352989375082 13202D0*x( 18) - 0.958821750950846804D0 + x(19) 

+ - 0. 8602462874967 7 1 303D0* x(20) 

y(4) = y(4) - 78.95852907402 15(X)4D0*x(21) + 

+ 43.922 1 7 34285092698DO*x(22) + 8.352989375082 13202D0*x(23) 

9 + - 1954.93624505O99399D0*x(24) - 

+ 2.434 1670262402 1802D0*x(25) - 2597.51251 185656298D0*x(26) 

, + + 37578.7624927093302D0 J » , x(27) - 

+ 4.17208401208873203D0*x(28) + 3.077920932(K)968897D9 ,, ‘x(29) 

+ - 41.0633443973 176799D0*x(35) + 

+ 3.573 16720948225297D-3*x(36) 

y(5) = -42.1666883238785903DO*x(5) + 1 6.3658705 1 303695D0*x(6) 

+ - 453. 129586 109553202D0 i, ‘x(7) + 35 1.46043821 35 13001D0*x(8) 

+ - 0.59975463706 1799697DO*x(9) - 

+ 4.63233083 144964902D-2*x( 10) - 1.5201 16146415901D-2*x(ll) 

+ - 6.863707281 8449770 lD-2*x( 12) + 

+ 0.48369850553 1 560 102D0*x( 1 3) - 

+ 0.969 166 155996266<)99D0*x( 14) + 

+ 0.519949869945893497D0*x(15) + 1. 2630282884686470 lD0*x( 16) 

+ - 1.74000709591828(M)lD0*x(17)- 

+ 1.1 89239579268 164D-2*x( 18) - 0.969166155996266099D0*x(19) 
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+ + 0.5 19949869945893497D0*x(20) 

y(5) = y(5) + 1.26302828846864701D0*x(21) - 
+ 1. 74000709591 82800 1 DO*x(22)- 1.1 89239579268 164D-2*x(23) 

+ + 900.378460475793005DO*x(24) + 

+ 0.643799490067852503D0*x(25) + 1238. 14470008882901 D0*x(26) 

+ + 29 1 . 10767237057 1 898DO*x(27) - 

+ 9.35961201 17822 1493D-2*x(28) + 3.81608823019032497D7*x(29) 

+ - 0.6521 1 6307 1 84986296DO*x(35) - 

+ 2.32409992594539401 D-5*x(3o) 

y(6) = 6384.77 199999999698D0*x(6) - 2628 1.1 30000000000 lD0*x(7) 

+ - 44252.8600000000997DO*x(8) 

y(7) = 8 1 .7875899999999092D0* x(6) - 1061.036D0*x(7) - 
+ 2228.2*D0*x(8) 

y(8) = -23 1 .7430000000004D0*x(6) - 3091.725999999995D0*x(7) + 

+ 61 62.58 100000000604D0*x(8) 

y(9) = 28.7296599999999702D0*x(6) + 534.1 0220000000 1204D0*x(7) 
+ - 5791.63599999999894D0*x(8) 

y( 10) = x(24) 
y(ll) = x(25) 
y( 12) = k { 21 ) 
y(13) = x(28) 
y( 14) = x(35) 

y( 1 5) = -2846.47376 1074609D0*x(5) - 2699.5767200 1826D0*x(6) + 

+ 6978.2941 1601534196D0*x(7) + 2. 14350293 155008902D5*x(8) - 

+ 22.32897639007 12601D0*x(9) - 3.3990629437 14242D0 + x( 10) - 

+ 9.13144768834393994D0*x(l 1) + 5.05173942454337099D0*x(12) 

+ + 75.7538 109829206405D0*x( 13) - 

+ 0.95882 1 7 50950846804D0* x( 1 4) - 

+ 0.86024628749677 1303D0*x( 15) - 78.95852907402 15004D0 + x( 16) 

+ + 43.922 1734285092698D0*x( 17) + 

+ 8.352989375082 1 3202D0*x( 1 8) - 0.958821750950846804D0 + x(19) 

+ - 0.86024628749677 1303D0*x(20) 

y( 15) = y(15) - 78.95852907402 15004D0*x(21) + 

+ 43.922 1 7 34285092698D0* x(22) + 8.352989375082 13202D0*x(23) 

+ - 1 954.93624505099399D0 * x(24) - 

+ 2.434 1 670262402 1 802D0*x(25) - 2597.51251 185656298D0*x(26) 

+ + 37578.7624927093302D0*x(27) - 

+ 4.17208401208873203D0*x(28) + 3.07792093200968897D9*x(29) 

+ - 41.0633443973 176799D0*x(35) + 

+ 3.573 16720948225297D-3*x(36) 

y(16) = -42.1666883238785903DO*x(5) + 16.3658705 1 303695D0*x(6) 
+ - 453. 129586 109553202D0*x(7) + 35 1.4604382 1351 300 lD0*x(8) 

+ - 0.59975463706 1799697D0*x(9) - 

+ 4.63233083 144964902D-2*x( 10)- 1.5201 16146415901D-2*x(l 1) 

+ - 6.863707281 8449770 lD-2*x( 12) + 

+ 0.48369850553 1560l02D0*x(13) - 

+ 0.969 166 l <; 5996266099D0*x( 14) + 

+ 0.519949869945893497D0*x(15) + 1 .2630282884686470 lD0*x( 16) 

+ - 1 .7400070959 1 82800 1 DO*x( 17) - 

+ 1. 189239579268 164D-2*x( 18) - 0.969166I55996266099D0*x(19) 

+ + 0.5 19949869945893497D0*x(20) 

y(16) = y( 16) + 1. 2630282884686470 lD0*x(21) - 
+ 1 . 7400070959 1828(X)lDO*x(22)- 1.1 89239579268 164D-2*x(23) 

+ + 900.378460475793005D0*x(24) + 
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+ 0.643799490067852503D0*x(25) + 1238.1447(XX)8882901DO*x(26) 

+ + 291. 10767237057 1898DO*x(27) - 

+ 9.35961201 17822 1493D-2*x(28) + 3.8l608823019032497D7*x(29) 

+ - 0.6521 16307 184986296D0+x(35) - 

+ 2.32409992594539401 D-5*x(36) 

c 1 - u bounded limit 

c — {nlarv/2.evel limiter.4} 

y(17) = min( 0.0174499999999999999D0, max( 

+ -0.0 1 7 4499999999999999 DO, y ( 1 ) ) ) 

c state-space system 

c -- {nlarw2.controller.l } 


c — ss ab - 

xd(l) = -27.7971 20754624 1899D0*x(l) - 25.4156325927014D0*x(2) 

+ + 0.966874297031 38 1502D0*x(3) - 

+ 8.71552238769890408D-2*x(4) - 9.51654087303377394D-3*yl4) 

¥ - 1.6803237 14059 14D-2*y(5) 

xd(2) = -6.87602622363934302D0+XH) - 
+ 7.367 1 1056293708599D0*x(2) - 8.27 1 35760534558596D-3*x(3) 

+ + 0.976 139521 572 193303D0+x(4) - 

+ 2.63798849248624(X)lD-3*y(4) + 1.501039866773363D-2*y(5) 

xd(3) = -52.3804248827 1 16199D0*x(l) - 
+ 35.5899410486724701D0*x(2) - 0.29(X)86695207493903D0*x(3) 

+ - 0.2 1 8 1482 1 3265887DO*x(4) - 7.09094969390322604D-3*y(4) 

+ - 9.45879838024 16<'603D-3*y(5) 

xd(4) = -4067.66699224439799D0*x( 1 ) - 
+ 507 1 .334 17368828702D0*x(2) - 5.043376280250499D0*x(3) - 

+ 17.025358 1968938801D0 + x(4) - 1.157405 14598378 101 D-2*y(4) 

+ + 1 .73620839739 127D-2*y(5) 

c state- space system 

c -- {nlarw2.arw2 moD.3} 
c — ss ab-- 

xd(5) = x(10) 
xd(6) = x( 1 1 ) 
xd(7) = x(12) 
xd(8) = x(13) 

xd(9) = -1 25.245923867737901 D0*x(5) - 
+ 45.03350413 19353904D0*x(6) - 1 17.872 130364775799D0*x(7) + 

+ 3906.6306 1647083703D0*x(8) - 1.26414087776252901D0*x(9) - 

+ 0. 1 27 3297 457 868026D0 * x( 1 0 ) - 0.200 1582942733044D0*x(l 1) + 

+ 0. 1463652541 766454DO*x( 12' + l.l 5969827063 154399D0*x( 13) 

+ - 0.969152580431533295D0*x(14) + 

+ 8.27755253907873397D-4*x(15) + 

+ 6.04388965594990202D-3*x(16) - 

+ 6.856145 16592095201D-3*x( 17) + 

+ 5.01807022310764295D-4*x(18) - 

+ 0.96915258043 1533295DO*x( 19) + 

+ 8.27755253907873397D-4*x(20) 

xd(9) = xd(9) + 6.04388965594990202D-3*x(2 1 ) - 
+ 6.856 145 1659209520 lD-3*x(22) + 

+ 5.01807022310764295D-4*x(23) + 1223.06162596226201D0*x(24) 

+ +0.899244 10982273 1702D0*x(25) + 

+ 1763.75969 19 1059301D0*x(26) + 698. 16 1055555356 lD0*x(27) - 

+ 0. 1264279353258893DO*x(28) + 6.7 1666950223660506D7*x(29) 

+ - 1. 867833656 l08772D0*x(35) + 7.6556232257 1364598D-5*x(36) 
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+ - 4.84533 1 1 555 1499498D-4+U 

xd(10) = -3 1.7973764346982 102D0*x(5) - 16.6021798016908D0*x(6) 

+ - 135.2985761432064D0*x(7) + 2147.94963514691301D0*x(8) - 

+ 1 .24869739759363099D0* x(9) - 0.2655459480335534D0*x(10) - 

+ 0. 153253792353896599 DO*x( 1 1) - 1.050858874417987D-2 + x(12) 

+ + 1 .278272593746799D0*x( 1 3) - 1.27530495555126799D-3*x(14) 

+ - 5.79 1 25 1 354 1 2297698DO*x( 1 5) + 

+ 2.918827546093392D-3*x(16) - 1.20106974790340501D-3*x(17) 

+ + 3.667313377607812D-3*x(18) - 

+ 1 .27530495555 126799D-3*x(19) - 5.79125135412297698D0*x(20) 

+ + 2.918827546093392D-3*x(21) 

xd(10) = xd(IO) - 1 .2010697479034050 1 D-3*x(22) + 

+ 3.6673 133776078 12D-3*x(23) + 4267.34662588543097D0*x(24) 

+ + 3.51 37358075365630 lDO*x(25) + 

+ 6539.00920066869003D0*x(26) + 779.8473335337912D0*x(27) - 

+ 0. 10065335997 10573D0*x(28) + 5.69431331591210403D7 + x(29) 

+ - 0.5308573 1374576990 lD0*x(35) + 

+ 6.99578 120780643495D-5*x(36) - 4.427709625 193978D-4*u 

xd(l 1) = -4035.453860038237D0*x(5) - 5039.804988 19288903D0 + x(6) 
+ - 103 15.3494347062699D0*x(7) + 2.4223448237 1904D5*x(8) - 

+ 27.9029350865074499D0*x(9) - 5.(X)483961999364602D0*x(10) 

+ - 1 6.9 1 99342752880302D0*x( 1 1) - 

+ 1.21 7457408 18824499D0*x( 12) + 75.01 44208438646292D0*x( 13) 

+ + 8.24300236275332708D-3*x( 14) + 

+ 1.93965721381322799D-3*x(15) - 152.11 3429 1940743D0 + x( 16) 

+ - 5.8038654 1425226303D-2*x( 17) + 

+ 5.25523764977082796D-2*x(18) + 

+ 8.24300236275332708D-3*x(19) + 

+ 1 .9396572 1381 322799D-3*x(20) 

xd( 11) = xd( 11)- 152.11 3429 1 940743D0*x(2 1 ) - 
+ 5.80386541425226303D-2*x(22) + 

+ 5.25523764977082796D-2*x(23) - 3324.495306375349D0*x(24) 

+ - 3.66 154656 177599203D0*x(25) - 

+ 3889.2879735487080 1 D0*x(26) + 70121.27667658 17 105D0*x(27) 

+ - 16.230958826885()602D0*x(28) + 

+ 8.2431213221 1215198D9*x(29) - 58.8147948818334498D0*x(35) 

+ + 6.03774882724283902D-3*x(36) - 

+ 3.8213600172423 1499D-2*u 

xd(12) = 387 1 .7533554938 1 299D0*x(5) + 1535.837069757428D0*x(6) 
+ - 34825.31 67988588 102D0*x(7) - 1 .03428353574 103699D5*x(8) 

+ + 26.295468997 1899098D0*x(9) + 1.2509891 7029 1227D0*x( 10) 

+ + 2.630979503694448D0*x( 1 1) - 18.3240096889994701D0*x(12) 

+ + 9.72678577610167805D0*x(13) - 

+ 4.201 120788 1 598 1 804D-3*x( 14) + 

+ 1.3530565 1065 17070 lD-3*x( 15) + 

+ 4.063760 124791 73 199D-3*x( 16) - 140.368755 107959299D0*x(17) 

+ + 0. 1 4 1 684 1 344 1 3 1 55902D0* x( 1 8) - 

+ 4.201 12078815981804D-3*x(19) + 

+ 1 .35305651065 1 7070 1 D-3*x(20) 

xd(12) = xd<12) + 4.063760 12479 173 199D-3*x(21) - 
+ 140.368755 107959299D0*x(22) + 0.141684 13441 3 155902D0*x(23) 

+ + 2449. 1 22483 1 5986699D1 )* x(24) + 

+ 2.20529981 83453 1003D0*x(25) + 1948.391 1057339500 lD0*x(26) 

+ + 20503.23546981 l2499D0*x(27) - 
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+ 

8.7 17721 34 128349308D0*x(28) + 3.55976045603601098D9*x(29) 


+ 

+ 56.3946 108862526 199D0+x(35) - 


+ 

3. 6337807 67 3383631 )2D-3 *x(36) + 2.29986124515086901D-2*u 
xd(13) = -463.806 1 800728 10101 D0*x(5) - 


+ 

349.7368901 3 1978697DO*x(6) - 1 195.977582526502D0*x(7) - 


+ 

40097.8732828805 196DO*x(8) - 1.46494257127194799D0*x(9) + 


+ 

0.1 83876792527 185799D0*x( 10) - 


+ 

0.5551 17409039766599D0*x(l 1) + 


+ 

0. 167 10 1 54453369 1 599D0*x( 1 2) - 29. 83280229208 1 8499D0*x( 1 3) 


+ 

- 3.8 1290566534986402D-3*x( 14) + 


+ 

5.78 1507063 19 158894D-4*x(l 5) - 0.1 0025484808 18097D0*x( 16) 

# 

+ 

+ 9.99420602 1040908D-2*x( 17) - 6.49269676792 184203D0*x(l 8) 


+ 

- 3.8 1290566534986402D-3*x( 19) + 


+ 

5.78 1507063 19 158894D-4*x(20) 

1 


xd(13) = xd(13) - 0. 10025484808 1 8097D0*x(2 1 ) + 


+ 

9.99420602 1040908D-2*x(22) - 6.49269676792184203D0*x(23) 


+ 

+ 1019.796982 139836D0*x(24) + 0.976905007059 1091 96D0*x(25) 


+ 

+ 2093.9417 1462465602D0*x(26) - 


+ 

5232.10908(K)6390696D0*x(27) - 1.29605906458 1286D0*x(28) + 


+ 

1. 285701 909877968D8*x(29) - 7.1 1996318 155942698D0*x(35) + 


+ 

6. 163729505 18372805D-4*x(36) - 3.90109462353402403D-3*u 
xd( 14) = 0. 1 2848 1 999999999999D0*x(9) - 


+ 

40.1736000(MMHKM)8()4D() + x(10) + 4.33304999999998597D0*x(l 1) 


+ 

- 98.3560999999999694D0*x( 1 2) + H)44.980000000003D0*x(13) 


+ 

- 84.254265530464 1794D0*x( 14) - 


+ 

146.22500(XHKKHX)399D0*x(25) + 158.846999999999799D0 + x(28) 


+ 

+ 8.527042966739 18397D-2*x(35) + 


+ 

4.06559376838870898D-2*x(36) - 0.25731 606 12904246D0*u 
xd(15) = -19.2853(XXXXXXXX)101D0*x(9) - 


+ 

5.5699 1999999999599D0*x( 10) + 10.0029999999999899D0*x(l 1) 


+ 

- 29. 13 120(XXXXXXX)4()2D0*x( 12) - 


+ 

92 . 8 560999999999694D0* x( 1 3) - 84.254265530464 1794D0*x( 15) 


+ 

- 1 1 1.353(XXXXXXXXX)99D0*x(25) + 


+ 

1 7 .2043999999999602D0* x(28) + 4.65810662194168197D-2*x(35) 


+ 

+ 2.2209304360874 1999D-2*x(36) - 


+ 

0. 1405652 1747388740 1 D0*u 



xd(16) = 6.8972 l(XXXKXXXXX)95DO*x(9) - 


+ 

6.32039(XXXXXXXX)295D()*x(10) - 5. 16849999999999499D0*x(l 1) 

J 

+ 

+ 2.0769999999999980 1 D0*x( 1 2) + 


+ 

205.770999999999699D0*x( 13) - 84.2542655304641794D0*x(16) 


+ 

+ 1 9 065 1 999999999999D0* x(25) - 

X 

+ 

59.3106999999999998D0*x(28) + 3.427498257283779D-4*x(35) 


+ 

+ 1.634190845563998D-4*x(36) - 1.034 298003521 5 19D-3*u 
xd(17) = -7.4448 19999999993DO*x(9) + 


+ 

6.89347(XXXXXXXX)798D0*x(10) + 8.567650(XXXXXX11509D0*x(l 1) 


+ 

+ 2.0852D0*x(12) - 145.588999999999899DO*x(13) - 


+ 

84.254265530464 1794D0*x( 17) - 15.7273(XXXXXXXX)lD0*x(25) - 


+ 

32.45039999999995( )4D0*x(2 8) + 3.19475760574534596D-2*x(35) 


+ 

+ 1.523222841035782D-2*x(36) - 


+ 

9.64065089263 157499D-2*u 



xd(18) = - 13.323899999999980 lD0*x(9) - 


+ 

18.16520(XXXXXXX)302D0*x(10) - 1 18.127000000000001D0*x(H) 


+ 

+ 1 34.6 1 4999999999799D0*x( 1 2) - 


+ 

949.754CXXXXXX)008(} 1 D0*x( 1 3) - 84.254265530464 1 794D0 + x( 18) 
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+ 

+ 3.1790500(XXXXXX)398DO*x(25) - 


+ 

397.355999999999803D0*x(28) + 0.27 8893228004259 lDO*x(35) 


+ 

+ 0. 1 3297 30 1 00 1 563080 1 DO* x(36) - 


+ 

0.84 160132921 2856797D0*u 
xd(19) = 82.2208000000000503D0*x(9) + 


+ 

38. 1446000000000804 DO*x( 10) - 9.01630999999997607D0*x(ll) 


+ 

+ 31.331 5999999999799 D0*x( 1 2) + 


+ 

656. 12 1999999999403D0*x( 1 3) - 168.6879863655722DO*x(19) + 


+ 

43 1 ,763000000(XX)801 DO*x(25) + 123.543000000000101D0*x(28) 


+ 

+ 0.7835 19507974698798D0*x(35) + 


+ 

0.373572883525770998D0+x(36) - 2.3643853387707 !699D0*u 



xd(20) = 2.01099(XXXXXXXX)699D0*x(9) + 

\ 

+ 

29.2762000000000202D()*x( 10) - 12.3716999999999799D0*x(ll) 


+ 

+ 20.6739(XXXXXXXXXX)2D0*x(12) + 


+ 

247.386999999999698D0*x(13) - 168.6879863655722DO*x(20) + 

1 

+ 

264.5491XXXXXX)00902D0*x(25) + 26.5588999999999902D0 + x(28) 


+ 

+ 1. 0063890730269250 lD-2*x(35) + 


+ 

4.79834470 1 47791 304D- 3 *x(36) - 3.03692702625 1846D-2*u 
xd(21) = 22.594(XXXXXXXXX)5D0*x(9) - 5.37708000000000697DO*x(10) 


+ 

+ 8.3899499999999989 lD()*x(ll) - 18.90300000000002D0 + x(12) 


+ 

+ 142.493(XXXXKXXX)4D0*x(13) - 168.6879863655722D0*x(21) 


+ 

- 22. 1969(XXXXXXXX)30 1 D0*x(25) + 


+ 

168.935999999999702D0*x(28) + 0.223846017964695498D0*x(35) 


+ 

+ 0.106727 147883002501D0*x(36) - 


+ 

0.675488277740527096D0*u 
xd(22) = -23. 13850001 XXXXX120 1 D0*x(9) + 


+ 

4.77663999999998601D0*x(10) - 12.28980000000001D0*x(ll) + 


+ 

13.8938999999999699D0*x(12) - 266.94 1999999999098D0*x( 13) 


+ 

- 1 68.6879863655722D0*x(22) + 19.15959999999995D0*x(25) - 


+ 

43.06591 XXXXXXXX)596DO*x(28) - 0.23 1364769970956501 D0*x(35) 


+ 

- 0.1 1031200038368 1901D0*x(36) + 


+ 

0.698 1772 176 1 824 l(X)3D0*u 
xd(23) = 1 50.079999999999899D0* x(9) - 


+ 

14.4340999999999999D0*x( 1 0) + 184.738999999999599D0*x(l 1) 


+ 

- 262.1530fXHXXXKX)198D0*x(12) + 


+ 

2774.5701XXHXHKX)698D0*x(13) - 168.6879863655722D0*x(23) + 


+ 

28.3 123tXXXXXXXX)5O2D0*x(25) + 935.655999999998997D0*x(28) 


+ 

+ 0.30235392999 1241402D0*x(35) + 


+ 

0.144158796714765301D0*x(36) - O.9123974475618084O6D0*u 
xd(24) = x(25) 

V 


xd(25) = -3.94784 179999999702D5*x(24) - 

t 

+ 

1 256.637001XXXXX ) 1 99D0*x(25) + 7.89568359999999404D6*x(26) 
xd(26) = -20.0D0*x(26) + y(17) 
xd(27) = x(28) 

xd(28) = -3.348D5*x(27) - 818.4(XXXXXXX)001498D0*x(28) + 

+ 

7.201548D+1 i*x(29) 
xd(29) = x(30) 

xd(30) = -2.151D6*x(29) - 540.099999999998502D0*x(30) + 


+ 

2.191D5*x(31) 
xd(31) = x(32) 

xd(32) = -2.191D5*x(31) - 185.3(XXXXXXX)0002D0*x(32) + 


+ 

3.742D6*x(33) 
xd(33) = x(34) 

xd(34) = -3.742D6*x(33) - 1446.5D0*x(34) + y(3) 
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t 


xd(35) = -0.3399999999999999 DO*x(35) - 
+ 0. 162 10799999999990 lDO*x(36)+ 1 ,02600000000000299DO*u 

xd(36) = -0.359999999999999397D0*x(36) + u 

c 

c 

return 

end 


3 


t 


« 
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Appendix D • MODEL. INP file for the ARW-2 model 


This is an example of the input file. MODEL.INP, that must be created by the user. The 
odd numbered lines are dummy variables used to name the fields and the even numbered 
lines contain the actual data. This is the input file used to generate the results shown in 
figure 5. 

model 

arw2 flexible nonlinear model, much .86, altitude 24k ft 
case 

maximize output #6 
nstates noutnoutmx 
36 17 6 

tmaximp deltat 

10.0 0.005 

sigma akmin akmax nkvals 

1530.0 10.0 15000.0 9 

iouttype impres iexcres iwave 
3 111 

standard filename (iouttype=l) 
arwl530.out 

matrixx filename (iouttype=2) 
arwl530.mat 


Where these input quantities are defined: 

model A character variable describing the model used in the analysis 

case A character variable describing the case 

nstates Number of states in the model 

nout Number of output quantities used in the model 

tmaximp The length of the simulation used in the impulse responses 

deltat The time step used in the simulations 

noutmx The number of the output quantity to be maximized 
sigma The gust intensity used in the analysis, units of velocity 

akmin The minimum k value to be used in the analysis 

akmax The maximum k value to be used in the analysis 

) nkvals The total number of k values to be used in the analysis 

1, akmin is used 

2, akmin and akmax are used 

f >3, k values distributed logrithimically between akmin and akmax 

iouttype Specifies the format for the output 

1, Normal ascii file 

2, MATRIXx readable ascii file 

3, Both 1 and 2 

impres If 1, write the impulse responses to the output file(s) 

iexcres If 1, write the excitation waveform responses to the output file(s) 

iwave If 1, write the excitation waveform(s) to the output file(s) 
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Appendix R - A sample listing of the MFB1DS program output using the input file 
shown in apppendix I) and the model listed in Appendix C 

nr\v2 flexible nonlinear model. Mach .86. Altitude 24k ft 
maximize output #6 

length of simulation for impulse responses = lO.OtXKH) 
number of time steps for inpulse responses = 2001 

length of simulation for excitation responses = 20.0000 

number of time steps for excitation responses = 4001 

gust intensity = 1530.00 

calculating impulse responses for each k value 
calculating normalized excitation waveforms 


1 

k = 

10.00000 sqrtt energy) of output 6 = 

568.177 

1 

k = 

24.9466 sqrl(cnergy) of output 6 = 

1417.29 

3 

k = 

62.2333 

sqrt(cncrgy) of output 6 = 

3536.37 

4 

k = 

155.25! 

sqrt(energy) of output 6 = 

8820.35 

5 

k = 

387.298 

sqrt(cnergy) of output 6 = 

22003.6 

6 

k = 

966. 1 77 

sqrt(cncrgy) of output 6 = 

56134.6 

7 

k = 

2110.28 

sqrt(encrgy) of output 6 = 

162952. 

8 

k = 

6012.84 

sqrl(energy) of output 6 = 

509979. 

9 

k = 

15000.0 

sqrttenergy ) of output 6 = 

1.4941 IH+06 

calculating 

maximized responses 


1 

k = 

10.00000 

maximum value of output 

6 = 2870(H). 

2 

k = 

24.9466 

maximum value of output 

6 = 286965. 

3 

k = 

62.2333 

maximum value of output 

6 = 286988. 

4 

k = 

155.251 

maximum value of output 

6 = 286997. 

5 

k = 

387.298 

maximum value of output 

6= 287025. 

6 

k = 

966. 1 77 

maximum value of output 

6 = 289885. 

7 

k = 

2410.28 

maximum value of output 

6 = 296994. 

8 

k = 

6012.84 

maximum value of output 

6 = 279944. 

9 

k = 

15000.0 

maximum value of output 

6 = 249730. 

saving data 
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